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Abstract 

A systematic construction of liigher order splines using two liierarchies 
of polynomials is presented. Explicit instructions on how to implement 
one of these hierarchies are given. The results are limited to interpolations 
on regular, rectangular grids, but an approach to other types of grids is 
also discussed. 
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1 Introduction 



The purpose of this work is to construct smooth interpolants for functions that 
are only known on the nodes of a regular rectangular grid. Cubic splines have re- 
cently been used in the context of particle or virtual particle tracking in complex 
fields, [1, 2, 3]. A more elaborate discussion on integrating particle trajectories 
in interpolated fields, and the advantages of using higher order splines, is to be 
published elsewere^. In this work just a review of the construction of the spline 
interpolations is presented, without rigurous proofs. 
Consider the functions of D variables 

/ : ^ M (1) 

and assume they can only be computed on the grid 

D 

G = n h,!^, (2) 

where < hj < 1 are the grid constants and hi, = {hz\z G Z}; a grid cell is 
defined as: 

D 

Basically, the values f{x) are only available if x e G — this can happen 
for a variety of reasons, the simplest being that they are measured from an 
experiment. In the following we will say that we are computing approximations 
of /. In a more rigurous setting we would say that we are building an array of 
spline functions that converge to a certain limit under certain conditions; when 
the function / is sufficiently nice (a term that still needs a clear definition), it 
will be equal to that limit. For instance, one of the properties of this limit is 
that any of its Taylor expansions converge everywhere (for example any function 
with a finite discrete Fourier representation). 

A polynomial spline is a function defined on M^, that is a polynomial on 
each cell and it has continuous derivatives everywhere up to a certain order. 
Note that the set of grid constants {/ij} is a given, and the limit spoken of 
before is the limit of very large orders of the polynomials entering the spline. 
This limit will always exist as long as / is well defined on the grid G; it is true 
that for certain functions the limits for different grids will not be equal, but 
these functions are irrelevant here. 

^C.C. Lalescu, B. Teaca, and D. Carati, "Implementation of high order sphnc interpolations 
for tracking test particles in discretized fields" , submitted for publication. 
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2 One dimensional case 



2.1 Hermite splines 

Assume that the correct values of a function f{x) are known on h'Z, and we 
are interested in finding an approximation for the interval [xq, xq + h] (the cell 
^{xq))- Without loss of generality, we express the variable x in h units, and the 
formulas will be deduced for the interval [0, 1] (the cell C(o))- 
On [0, 1] we construct the n-th order polynomial 

n 

.W(x)=^a("'x'= (4) 

For Hermite spline interpolation, it is imposed that, on the enclosing grid nodes 
s^"-* coincides with the original function and the derivatives of s^"^ up the order 
m = (n — l)/2 coincide with the derivatives of the original function: 

^(z) = /W(x)l l^o;7^ (5) 

ixe{o,i} 



dx' 

By solving the linear system of equations (5) the coefficients aj,"'' depending 



where we called the ^-th order derivative = -^f 

By solving the linear system of equatio 
on /(O), /(I), /'(O) . . . can be easily foimd 

m 1 

EE^i'?^/^'no- (6) 



(") _ 

1=0 i=0 



Here we will discuss a method that avoids the computation of these coefficients. 
This leads to rewriting the expression of the spline as: 



ml / n 

= EE/''HqE^, 



kli 



S 



1=0 i=0 \k=0 / 

m 

^"^(-) = EE/'H*K'"'(-) (7) 

;=o i=o,i 



where the a polynomials can be found for specific values of n by symbolic 
computation; they can be defined equivalently as the solutions of the following 
system of equations: 

^ai"^'°)(x)=4,5,,, 1 = 0;^ (8) 

j,tG{0,l},x=] 

The set is a hierarchy of spline approximations (Hermite splines) , and 

it is based on the spline polynomials of the first kind a["'^^ ■ These polynomials 
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have the following explicit expressions: 



; m — l , , 

4"'"(-) = y(i--r+^E^^^-' (9) 

k=0 

= ^-^.--^ "f (^(1 - (10) 

fc=0 

(with aj"''^(x) — (— l)'ag"''''(l — x)); see appendix B for a proof. 

The distance 118^"+^^ — s^"^|| will obvioulsy go to for large n, because the 
splines coincide with the Taylor approximation around the two nodes up to the 
order m; note that this distance can be defined as the maximum of the distance 
|g(ri+i)^2;) — s^")(x)|. And there is a class of functions / that are limits of such 
approximations (all piecewise polynomials for instance). 



2.2 Grid splines 

If centered differences are used to compute the derivatives, the formula can be 
further adapted. The result is called here a grid spline (it should be a type of b- 
spline), as it is found solely from the values of the function on the grid. Centered 
differences are used because generally for D > 1 dimensions — when discussing 
practical implementations of higher order splines — the memory cost of keeping 
all the derivatives necessary is prohibitive (they arc as many as the coefficients 
aj,"'). A finite difference (for this rcscalcd variable) is just a linear combination 
of the numbers f{z) with z G Z. Also, it is crucial that the use of centered 
differences imposes the continuity of the derivatives when the polynomials are 
put together. 

To be more specific, rewrite the Hermite spline as 

^^"H-)-EE^i'"^'(-)- (11) 

/=0 2=0,1 

The smoothness of the interpolant is only determined by the fact that the coeffi- 
cients of the Taylor expansions used, ■* , are determined by the grid node i (they 
are the same whether the node is approached from the left or from the right, or 
more simply they are invariant to the cell). In practice the centered differences 
are the coefficients of the Lagrange interpolation polynomial J2^=o ^(f g)-^'^' with 
the coefficients determined from the equations 



23 



.k=0 



_ (12) 

-a, a 



Because the information contained in this Taylor expansion must be invariant 
to the cell, it must be obtained from a set of grid nodes that is symmetrical to 
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the current grid node i (thus the equations are solved at a: = —g, g, so a unique 
solution is obtained for a polynomial of order 2g). 

It is then obvious that a centered difference f<i'^> « can be written as: 

/<^''>(^)^^«,= ^ c^'V(* + fc) (13) 

k=-g 

where the coefHcients c^"^''^ are numbers that only depend on g. This means 
that the final formula of the spline will use the values of the function on the 
q ^ 2g + 2 nodes —g, . . . , 0, 1, . . . , g + 1: 

1+5 

where the spline polynomials of the second kind /3|"'*^ (x) can be found for a 
fixed n and q by symbolic computation. 
We give the pf'^^'-s as an example: 

p'-^C\x) = \{x~lfx{2x+l) (15) 

f}f-'^){x) = ~]^{x-l) {Qx^ ~Qx^ + 2x + 2) (16) 

Pi'^^ {x) = \^ (6 •t'' - 15 a;^ + 9 + x + l) (17) 

f}^^''\x) = ~\[x-l)xH'2x-i) (18) 

In terms of accuracy, using the distances 1 1/<«''> - /^'^ 1 1 one should be able 
to find the distance ||s("'«) - s(")||. 

As a sidenote, it is quite clear that the same approach can be used for 
irregular grids (the Taylor expansions can still be approximated using adjacent 
nodes). However, some of the simplifications that can be made for regular grids 
will no longer be possible. 



3 D dimensions 

Once the spline polynomials of the first and the second kind are available, the 
spline interpolation formula can be directly generalized to a scalar function of 
D variables (expressed in hj units, and shifted to [0, 1]^): 



s(-\x,,x,,...^xn)^ /('-•^'-)(*i,...,zz,)n"!:'^^( 



D 

li,...,ln=Oii,...,io=0,l k=l 

(19) 

i+g D 
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A rather interesting observation is that for the Hermite splines exactly 2^ (m + 
1)^ input values /^'i^ - ''"^ (ii, . . . ,io) are required for a given n = 2to+1, while 
for grid splines input values are required for a given q (and m < 2g, thus 
n < 2(7 — 3). This means that grid splines generally achieve a given degree 
of smoothness from less information than a Hermite spline — and the price is 
probably a much larger error. Also, the number of input values is the number 
of terms in the sum to be computed, thus grid splines will be faster to compute 
than the corresponding Hermite splines (up to a maximum q that will depend 
on the order). As a final note, the mixed derivatives of these fields are also 
smooth: 




(21) 



(whenever m > max{lj}). 

4 Implementation of the grid splines 

In practice the Hermite splines will probably not be very useful, as they require 
too much memory. Other than that, their implementation should be similar 
to that of the grid splines. Note that we mention "parallelized codes" in the 
following; this refers specifically to cases of computer programs that run on 
several processors at once, with the memory divided between them, and these 
programs work with physical fields, each processor keeping a slice of these fields 
in its memory. 



4.1 Full expressions of the ID grid splines 

You will need to use a computer algebra system. For q = 4:, define the following 
functions: 



(22) 

(23) 

(24) 

(25) 

(26) 
(27) 





- f" 
" ^(0,1) 




2 






— y^" 




-V,i)^ 






~ *0 


, (3,4) 


, (3.4) 2 


Vsf'^x^ 




- ,(5'4) _ 
~ *0 


1 (5,4) 


, (5.4) 2 
{■ S2 X - 


K 45^4)^3 ^^(5,4)^4 ^^(5,4)^5 



Afterwards, solve the following systems of equations (symbolically): 

fi(o,i)(0) = /(0), t(o,i)(-l) =./(-!), t(o,i)(l) =/(!), 
U(i,i)(0) = /(1), t(i.i)(-l)=/(0), t(i.i)(l)=/(2). 



,(3,4) 



(0) 



" d 

dx 



= (3,4) 



{x) 



x=0 



'(0,1) 



'•(0,1) 



= (3,4; 



(1) 



' d_ 

dx 



= (3,4) 



{x) 



x = l 



"(1,1)' 



'(1,1)' 
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= (5,4) 



(0) 



ax 



x=Q 



x=0 



"(0,1) 



"-(0,1) 



"-(0,1) 



= (5,4) 



(1) 



ax 



ni,l)' 



ni,l)' 



(28) 



You can either solve them one at a time (and afterwards replace the solution 
of (26) into the expressions of the splines), or you can solve them all at once. 
For larger numbers of grid points q, the process is identical. 



4.2 Spline polynomials 

Note that the same systems of equations need to be solved for the Hermite 
splines, just that the centered differences f^^- should be replaced with the /^'^(i), 
and the spline polynomials of the first kind can be found as: 



a. 



("■')(^) = ^^^^(n)(^) (29) 



d(/(0(z)) 



For the spline polynomials of the second kind, the values of the function come 
into play (and the following expression makes sense after the centered differences 
have been introduced into the expressions of the splines): 

These two equations arc the simplest way to find the polynomials. Considering 
that the full expressions of the splines are linear in the t|j ^yS which are linear in 
the f{i)-s, they are equivalent to rearanging the terms in the sum and extracting 
the "coefficients" of the /(i)-s, as in the definitions. 



4.3 Implementation 

After finding the spline polynomials of the second kind, they should be put into 
their Horner forms (for fast computation). The following steps depend on the 
programmer's taste and abilities mostly, but we recommend the implementa- 
tion of the subroutines that follow. They are easy to adapt for the case of a 
parallelized code, and this implementation proved to be quite efficient and easy 
to work with (easy to expand for more cases, explain to other users, check for 
errors when necessary). Note that the same structure can be used for interpo- 
lating derivatives of the field if necessary, just that subroutines for computing 
the derivatives of the spline polynomials have to be added. 

Recommended structure of the interpolation code (subroutines): 

1. spline polynomials (ti, q) 
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• input: fraction ^ 

• output: 7i = (array of dimension q) 

2. grid coordinates 

• input: "normal" point coordinates {xi, . . . , xd) 

• algorithm: compute each Xj = [xj/hjl and each xi = xj/hj — Xj. 

• output: 

— the set of integers (ii, . . . , xd) 

— the set of fractions {xi, . . . ,xd) 

3. spline formula 

• input: 

— the set of fractions {xi, . . . ,xd) 

— the type of spline (ti, q) 

— a pointer (or similar notion) / to an array containing the informa- 
tion about the local field (the values of the field on the nodes of 
the cell C'(ii....,4o) and the necessary neighbouring cells), shifted 
such that /(0, 0, . . . , 0) = /(xi, . . . ,xd)- 

• algorithm: compute the polynomials /J^^"'*'' (ij) in the array 7y (by 
calling the spline polynomials subroutine), then compute the sum 

l+g D 

f= E (/(H,...,*D)n7y); (31) 

ii,...,ir>=0-g j = l 

note that testing shows it is more efficient to introduce as little do 
while loops as possible — for our 3D implementation, for q = 4: the 
sum is written in full in the source code, as introducing do while 
type loops slows down the code considerably. For higher values of q 
we just have one do while loop for one of the variables. 

• output: the approximation /. 

4. wrapper 

• input: 

— "normal" point coordinates {xi, . . . ,xd) 

— a pointer / to the array containing the information about the 
entire field 

— the type of spline (ri, q) 

• algorithm: put the local field in the array / from / and then compute 
the approximation / (using the above subroutines) 

• output: the approximation / 
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The wrapper is very useful in the case of a parallelized code. All the op- 
erations related to bringing together information spread on possibly several 
processors can be placed inside the wrapper, allowing for easy debugging and 
maintenance of the code. 

As an example of a parallelized version, in our implementation a 3D field is 
divided along the z coordinate between processors. The field is periodic in all 
directions, and obtaining / implies a little care in regards to the z coordinate, 
but it basically just requires a normal use of the MODULO operator. We impose 
that each processor has at least q nodes on the z direction in its memory, so 
that the formula contains at most information from two processors ("low" and 
"up"). We then compute (31) on each processor, only for the nodes that are 
in "its domain", and we then sum the two resulting values / = jiow + ]up\ the 
amount of information passed between processors is thus kept to a minimum. 



2g D 

ki 



A Appendix: Uniqueness 

Construction for general case: 

t(i,5)(u)= ^(i,.)!!"^ (32) 

ki,...,ko=0 i=l 
n D 

fel fer>=0 j = l 

with ii, . . . , G {0, 1}. 

The Taylor expansions (the centered differences) are given by the system of 
equations 

t(i,<,)(v) = /(i + v) (34) 
where ui, . . . , = — g, g. This system of equations has a unique solution t^j ss 
(there are as many t^^ as there are /(i + v), and this is in fact a simple 
Lagrange interpolation). 

The splines are given by the system of equations 



n(i-) I 



<(...) (35) 



which also has a unique solution (and q = 2g + 2), so formula (20) must give 
this unique solution. 

A generalization of this method would be to go to other types of grids. For 
instance one could imagine the case where in three dimensions there is a grid 
made up of equilateral triangles in the (x, y) plane, and squares in the other 
two planes. What changes is the fact that the sums are a bit more complicated. 
The crucial property that has to be preserved to have a smooth approximation 
is that the Taylor expansion must be the same no matter from which cell we 
approach a given node. 
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B Appendix: Spline polynomials 

For a rigurous analysis of the errors of these approximations, the distances 
!!/<«■'> ||s("'«) -s(")|| and ||s("^-/|| should be computed. After using 
the method presented above to find the a polynomials up to ti = 19, it can be 
seen that all these spline polynomials of the first kind have a simple form: 

«r(.-) = fr(i--r"^E^^-'= (36) 

fe=0 



-^^E^(i-)^ (37) 



l\ ^ m\k\ 

k=0 

(with aj"''^(x) = (-1)'4"''^(1 - x)). 

It would be useful to have an explicit formula for the a polynomials in the 
general it would allow for a more straightforward treatment of the error 

_ which is close to ||s(") - /||. 

Obviously, if this form is proven to apply for ag"''\ the form for a^"'''' must 
also be correct. 

First, notice that for x very close to 1 {x = 1 — y, with y small), the Taylor 
expansion of begins at y™'*'^, and this is half of the proof. To continue, 

note rewrite the polynomial as 

J^'i)(r)- ^ y^y^ + fc)! + l)!(-l)% ;+fc+,- 

"° ^^^-^Z^Z. fc! j!(^ + i__^)!^ ^"^^^ 

fc=o i=o ^ 

For this proof the coefficient of x'" in ag"'''', for < < ^ is needed. Obviously, 
for Iq < I this coefficient is 0: 

l + k + j = lo (39) 

k + j = lo-l (40) 

but /c + j > (41) 

so - ^ > (42) 

For Iq = I, the coefficient is given by the term with k + j = 0: 

1 (m + O)! (to+ 1)!(-1)0 _ 1 
n 77i!0! 0!(m + l-0)! ~ U 

which is what is needed. For Iq — I ^ Al > 1 we have: 

/ , An l^'"^'(m + fc)!(m + l)!(-l)i 
/ ■'^ m'.kl i m + l — 7) 



(43) 
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(with the Kronecker S) . This formula simplifies to 



lg(m + fc)! (m+ l)!(-l)^'-'= 

k=0 

And in fact this sum is, from [4]: 



5(m, I, AO = ^ ^ (A/~fc)!(™ + l-A/ + fc)! ^^^^ 



(-l)-^'sin(^A/) 

^^"^'^'^') = — n^A] — 

which is for integer values of A/ (note that A/ < m for the sum to make 
sense) . 
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